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We investigate dynamics of deformable self-propelled particles with a repulsive interaction whose mag- 
nitude depends on the relative direction of elongation of a pair of particles. A collective motion of the 
particles appears in two dimensions. However this ordered state becomes unstable when the particle den- 
sity exceeds a certain critical threshold and the dynamics becomes disorder. We show by a mean field 
analysis that this novel transition characteristic to deformability occurs due to a saddle-node bifurcation. 
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Dynamics of self-propelled objects have been studied for 
these more than one decade. Various collective motions have 
been found for interacting particles' where the hydrody- 
namic formulation has been developed."*'^* See ref.^' for a re- 
view of experimental studies. Effects of noises have also been 
investigated in various model systems to find a rich variety 
of collective dynamics.^ However, almost all of the previ- 
ous studies assume that objects are rigid despite the fact that 
coupling between shape deformation and migration in self- 
propulsion is common in biological microorganisms'^''^* . In 
the present paper, we are concerned with the deformation ef- 
fect on the collective dynamics of self-propelled particles. 

Quite recently, we have introduced the model system where 
the objects are deformable.'"'* We have shown that deforma- 
tion produces a lot of interesting dynamics of individual par- 
ticles such as rotating motion, zig-zag motion, chaotic motion 
and helical motion, '^-i?) These results clearly indicate that de- 
formation is relevant to the dynamics. In the present paper, 
we study dynamics of assembly of deformable self-propelled 
particles where there is a repulsive short range interaction be- 
tween a pair of particles such that the magnitude depends on 
the relative direction of elongation. Our main concern is sta- 
bility of collective dynamics due to deformation of individual 
particles. 

We consider interacting particles in two dimensions. The 
time-evolution of the i-th particle obeys the following set of 
equations for the position of the center of mass the veloc- 



ity v^^ and the deformation tensor S 
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where y > 0, k > 0, a and b are constants. This set of equa- 
tions is an extension of a single particle case.'"'' The tensor 
S is represented in terms of the unit normal n parallel to the 
long axis of an elongated particle as 
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with i, the magnitude. The last term in eq. (2) stands for the 
interaction between particles and is given by 



7=1 



(5) 



where is a positive constant, the total number of parti- 
cles and Fjj - -dU(rij)/drij with r,y - ri - rj. The form 
of the potential Uirij) will be specified later The other fac- 
tor Qij expresses the dependence of the relative angle of a 
pair of elongated particles such that when the interacting two 
particles are parallel, the repulsive interaction is weaker It is 
defined by 
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where Q is a positive parameter 

We have introduced the force (5) to take into considera- 
tion, approximately, of the excluded volume effect of elon- 
gated particles. Since the force depends on the variable 5^^, 
it is possible generally that there is a counter term in eq. (3). 
However, for simplicity, such a contribution is not considered 
in the present study. 

Equation (3) implies that if the velocity is zero, v[j' = 0, 
the particle does not deform and takes a circular shape in two 
dimensions. One of the characteristic features of the set of 
equations (1), (2) and (3) is that, when the interaction term 
is absent, these equations do not contain the particle size ex- 
plicitly as a parameter. Therefore the interaction range in the 
potential U is the only fundamental length scale of the system. 

We introduce an orientational order parameter defined by 
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where 9j is defined through the relation n*'' - (cos 6*, , sin ft). 
In order to quantify the degree of coherent motion, we may 
define an alternative order parameter'' 
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where we have put v*'' - |v'''| (cos 0,, sin 0,). In the following, 
we shall use <t> but we have verified numerically that the dif- 
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Fig. 1. (Color online) Snapshots of the time-evolution of particles at (a) 
t = 256 and (b) t = 512 for p = 0.075 and N = 8192 starting with a 
random initial configuration. These figures display only the 1/4x1/4 area 
of the system. The colors distinguish the direction of elongation and the 
arrows indicate the direction of migration. 
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Fig. 2. (Color online) Growth of order for (a) p = 0.0032 and (b) p = 
0.0852, and for N = 512. The curve given by eq. (10) is shown by the 
broken line in Fig. (a). 



ference between (S> and O,, is negligibly small in the plots in 
Figs. 2 and 3 below. 

We shall solve eqs. (1), (2) and (3) numerically. Hereafter 
we show the results for the potential 



U(rij) = exp 



(9) 



The parameters are chosen as a = 1 and b - 0.5 so that each 
particle tends to elongate to the direction parallel to the veloc- 
ity. ^'^^ We fix the parameters y - k - \ in which the particle 
undergoes straight motion.'"** Unless stated explicitly, other 
parameters are chosen as cr = \, K - 5 and Q - 50. The set 
of equations is solved by the Euler explicit method with the 
time increment 6t - 1/64. Particles are confined in a square 
area with size L and the periodic boundary conditions are im- 
posed. In most of the simulations, we choose N =512, but the 
cases =1024, 2048 and 8192 are also examined to check 
the number dependence. The density p = N/L^ is varied by 
changing the system size L. 

Appearance of the collective motion of self-propelled par- 
ticles is shown in Fig. 1 starting from a disordered initial con- 
dition. The particles constitute locally an approximate hexag- 
onal lattice such that one of the fundamental reciprocal lat- 
tice vectors is parallel to the propagating direction. Another 
orientation of the lattice such that one of the primitive lattice 
vectors is parallel to the propagating velocity has not been ob- 
served. We note that the propagating lattice is not completely 
stationary. Each particle exhibits a small oscillation (or fluctu- 
ation) in its center of mass perpendicularly to the propagating 
direction. 

Figure 2 shows the growth of the order parameter O for 
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Fig. 3. (Color online) Hysteresis loop for (a) N = 512, (b) N = 2048 and 
(c)A' = 8192 



the dilute regime p - 0.0032 and the intermediate regime 
p = 0.0852 obtained for N - 512. Initially the velocity and 
the position of the individual particles are random. When the 
density is low, the order parameter starts to grow at around 
t = 3000 and almost saturates at about t - 7000 as in Fig. 
2(a). When the density is high, the transition to the ordered 
state is more abrupt. Figure 2(b) indicates that there is an in- 
cubation time and the order parameter starts to grow at around 
t = 1500 and completes within a narrow interval. 

It seems that the collective motion appears asymptotically 
in time even for a dilute limit since once an ordered motion 
appears after collisions, there is no mechanism to destroy it. If 
one adds noise terms in the time-evolution equations, a sharp 
transition from the disordered state in the dilute regime to 
the ordered state is expected to occur However, we do not 
carry out such numerical simulations because there have been 
many studies for rigid particles with noises^""' and, further- 
more, the deformability is expected to be irrelevant in the di- 
lute regime. 

The time-dependence of the order parameter in Fig. 2(a) 
can be fitted by the following function 

(10) 

where a = 5.02 x lO"",/? = 3, y = 0.955 and fo = 2.25 x 10\ 
If fo is put to be zero, no satisfactory fit is possible. However, 
the abrupt change shown in Fig. 2(b) cannot be represented by 
this type of function. This implies that the transition kinetics 
is quite different from the ordinary nucleation and growth in 
first order phase transitions in thermal equilibrium. In fact, we 
have seen that once an ordered region appear in the present 
self-propelled system, it spreads over the system almost in- 
stantaneously when the density is moderately high. 

The ordered state becomes unstable for large values of den- 
sity. We have investigated the disordering process by increas- 
ing the density slowly. The compression and expansion are 
made at a constant rate 



0(0 = r[i 



where 



dp/dt - ±A, 



/I = (pi -po)2 



(11) 



(12) 



with typical values po = 0.07 and pi = 0.14. The expo- 
nent « is varied. Figures 3 display the order parameter O as 
a function of p for n = 14. It is evident that the ordered 
state is spontaneously broken when the density exceeds a cer- 
tain threshold value and there is a hysteresis in the transition. 
Figure 3 shows the A^-dependence of the hysteresis behav- 
ior The threshold pdo for the transition from the disordered 
state to the ordered state by decreasing the density is less 
sensitive to the system size. We note that the threshold for 
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Fig. 4. (Color online) The oder-disorder transition density p^j as a function 
of Q for N = 512. The dotted Une is a guide to the eye. The volume fraction 
of the total particles defined by ify^i = ncr-p^^ is also indicated as the scale 
of the vertical axis. 




1.0- 
0.8- 
.0.5- 
'0.4- 
0.2- 
0.0- 



10" 



10^ 



10^ 



10" 



Fig. 5. (Color online) Diffusive property of the disordered state for p = 
0.14 and W = 512. (a) Self-intermediate scattering function for wave num- 
bers kcr = 1/16, k<T = 1/8, k<T = 1/4, kcr = 1/2 and kcr = I from top to 
bottom. The dotted curves show the function given by eq. (15). (b) Mean 
square displacement averaged over 16 independent runs. 



the ordered-to-disordered state by increasing the density oc- 
curs at Pod ~ 0.132 for A? = 512 = 2^ pod ~ 0.125 for 
N = 2048 = 2" and Pod ~ 0.115 for = 8192 = 2'^ This 
indicates that the hysteresis region pod - Pdo has a tendency to 
decrease Hnearly with In A^. However, more accurate analysis 
is necessary to obtain any definite conclusion, which is left for 
a future study because of time consuming of numerical simu- 
lations for larger system. We have also carried out numerical 
simulations for = 5 12 for slower change of the density with 
n - lb and « = 18. The threshold density pod is approximately 
given by pod = 0.132 for n - 14, pod - 0.131 for « - 16 and 
Pod - 0.128 for « — 18. On the other hand, the threshold pdo 
is given by pdo = 0.0840 for n = 14, pdo = 0.0891 for n = 16 
and Pdo = 0.0902 for n = 18. Therefore, as expected, the 
hysteresis region tends to shrink for slower density change. 

The transition between the ordered and disordered states in 
the high density regime has not been observed in any rigid 
self-propelled systems and is characteristic to the deformable 
particles. The ordered state exists irrespective of the value of 
the anisotropic parameter Q in eq. (6). However, the stability 
limit of the ordered state depends on it. The critical density 
decreases as the value of Q is increased as shown in Fig. 4. 
When the anisotropy is absent, Q - 0, the transition to the 
disordered state does not occur. 

We have explored the properties in the disordered state at 
the high density regime by evaluating numerically the self- 
intermediate scattering function define by 

Fs{kj)^Uj^e->'^<^'('^-'-A. (13) 

The mean square displacement of the individual particle can 
be obtained from Fs{k, f) as 

C{t) = ^^^tJ^I, = {m - r(0)|2) = ADt, (14) 

where D is the diffusion constant in two dimensions. 

The intermediate scattering function is shown as a function 
of time in Fig. 5(a) for several values of A: = |^|. It is clear 
that there is no multiple relaxation which is typical in glass 
dynamics.'**' Actually, as shown in Fig. 5(a), the scattering 
function can be fitted by 

F,(^,f) = exp[-(fM)^'], (15) 

where ^ 23.46 x k^^ with v = 1.85. The exponent j3k in- 



creases from j6jt K 1.09to/Si K 1.3 by increasing A:cr from 1/16 
to 1 . This indicates that the usual diffusion with the wavenum- 
ber dependence fikV ~ 2 occurs for small values of k whereas 
a ballistic motion is dominant in the short length scale and 
clearly excludes the stretched exponential with yS^ < 1 . It 
is also verified (not shown) that the simple exponential decay 
independent of the values of k is less satisfactory to fit the nu- 
merical data. The mean square displacement is displayed in 
Fig. 5(b). This has been obtained by the average of 16 inde- 
pendent runs. It is approximately proportional to time asymp- 
totically with the magnitude of D a; 0.2. It is evident that there 
is no regime where C{t) has a plateau due to a cage effect. 
This also implies that the disordered state is different from 
the usual glass state'^' . 

We have shown that the coherent ordered state of the de- 
formable self-propelled particles becomes unstable for suffi- 
ciently large density and a disordered state appears. We have 
verified that the same transition occurs not only for the soft 
potential given by eq. (9) but also for 

U = exp(-r,ycr). (16) 

For example, the order-disorder transition occurs at around 
p = 0.02 for decreasing the density and p = 0.06 for increas- 
ing density for A^ = 512 keeping all the parameters the same 
values as used above. 

Here we describe a mean field analysis to clarify the order- 
disorder transition in the high density regime. The force term 
in eq. (2) can be written as 

N 

vC) . /"■) = /:v<'> ^ v^'* ■ FijQij = Kv^^N (v"' ■ FijQi) , (17) 

where the unit vector v^'* = v*'V |v*''| is regarded as a stochas- 
tic variable. The average is defined such that 

First, we make a decoupling approximation such that Qij in 
eq. (17) is replaced by the average {Qij), 

v«./') =/:v»A^(v*')-F,;)(e,v). (19) 

By a dimensional analysis, we note that A^ ^v''' ■ Fij^ might 
have a factor crp and have verified numerically that this quan- 
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Fig. 6. Correlation as a function of p for A' = 512. The mean values are 
indicated by the black dots whereas the standard deviations are shown by 
the error bars. This plot is obtained by decreasing the density slowly. The 
transition occurs at p = 0.082. The inset shows that^ is also negative for 
the ordered state. 

tity is indeed proportional to the density p so that we put 

A^(v«-f,y)=pa-;t. (20) 
The unknown quantity x is to be evaluated numerically. An- 
other decoupling approximation {(Sai3)~') - (s aj}) is also 
employed. Secondly, we eliminate S ap by putting dS apldt - 
in eq. (3). By using these approximations, eq. (2) becomes 

^ yw - gw- + X + Q(w - {v-'j)\ (21) 

where w = > with the superscript (;) omitted and 
g= 1+ (ab)/(2K), and Q = Qb^/{4K^) and x = Kpxo: Putting 
^v~^ = (v)^ + {(6\')~'^ with 6\' the fluctuating part, we iden- 
tify (v)~ with the stationary stable solution of eq. (21). We 
have verified numerically that;(f is negative as shown in Fig. 
6. Therefore, we note that w = is always a stable station- 
ary solution and that there is a pair of finite solutions, one 
of which is unstable and the other larger solution is stable. 
For sufficiently large values of \x\, this pair of solutions dis- 
appear. Therefore there is a saddle-node bifurcation at some 
value of p s Pod beyond which only the solution w = exists 
and hence this state corresponds to the disordered state. In the 
disordered branch, we should put (v)~ = in eq. (21) whereas 
the stable finite stationary solution should be equated to (v) 
in the ordered branch. The fact that < is essential for the 
saddle-node bifurcation. However, a quantitative comparison 
is impossible with the present crude approximation. For ex- 
ample, we have to require = -0.57 (x = -0.73) to realize 
Pod = 0.12 (pdo - 0.082), which are substantially different 
from the values in Fig. 6. Note also that the coefficient of the 
w^'- term which is dominant for large values of w in eq. (21) 
has a factor pQ. This implies that, if x and {{6v)~'^ are not 
strongly dependent on Q, the threshold pod decreases as Q is 
increased, which is consistent qualitatively with Fig. 4. 

We mention that the above mean field theory does not pre- 
dict the intrinsic hysteresis behavior of the transition in the 
slow limit of the density change and the N ^ oo limit. 
We need to take account of the stochasticity which is self- 
generated from the nonlinear interaction among the self- 
propelled particles. 

We emphasize that this kind of order-disorder transition has 
not been observed, to our knowledge, in any rigid rod-shaped 
particles and hence it is a characteristic feature of deformable 
soft particles. It is mentioned here that nonequilibrium phase 



transitions due to a saddle-node bifurcation have been studied 
theoretically in a single-variable Langevin system.-*" 

Finally we discuss some apparently related phenomena. 
Colloidal particles (although not self-propelled) exhibit glass 
behavior for high density. In this glass dynamics, motion of 
each particle is frozen or trapped in a cage and does not un- 
dergo a usual Brownian motion. This is contrast to the present 
result in Fig. 5 where particles have a finite diffusion con- 
stant though the origin of the motion is not thermal but comes 
from the self-propelled energy. A jamming transition occurs 
in granular materials under shear flow. The diffusion constant 
is finite in the jammed state but it is defined perpendicularly 
to the flow direction-'' and inherently anisotropic since the 
system is subjected to shear flow. The third related system is a 
traffic flow which also exhibits a jamming state when particles 
are dense. However, the flow is generally directed in one di- 
mension^^' or anisotropically one-dimensional even in the 
two-dimensional space^"*' and therefore the jammed state has 
no direct analogue with the present isotropic disordered state. 
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